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Abstract 

The pair distribution function of the electron gas is calculated using a parameterized gener- 
alization of quantum hypernetted chain approximation with the parameters being obtained by 
optimizing the system energy with a genetic algorithm. The functions so obtained are compared 
with Monte Carlo simulations performed by other authors in its variational and diffusion versions 
showing a very good agreement especially with the diffusion Monte Carlo results. 
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I. INTRODUCTION 



The Sommerfeld electron gas model[T] has proved to be very useful in describing many of 
the electronic properties of metallic solids. It represents the conduction electrons as a zero 
temperature ensemble of point charged fermions moving against a continuous neutralizing 
background that plays the role of the ionic lattice. In the simplest version of the model, 
fermions are considered spinless and the background is just characterized by a dielectric 
constant. If we have N fermions of mass m each one carrying a charge e, then the system 
Hamiltonian reads 

^ ' i=\ i<j 

Here v {vij) is the pair potential given by 

<r^j) = ^^ (2) 
where rjj = |rj — rj| with r,; the position of the ith-particle and e denoting the dielectric 
constant. 

The equilibrium behavior of this system has been widely studied through quantal Monte 
Carlo simulations [2], [3] and also from a variety of many-body theories [1]-[H]. Many of them 
center on the pair distribution functions (PDF) [9]. If the we denote ip {ri ,r2 , ■ ■ ■ , tat ) the 
A^-body wave function then the pair distribution function is given 



p (ri , r2 ) = p (ri ) p (rg ) fi- (ri , 

[... [dr 

= N{N-1)^- ^- (3) 
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where the integrations are over the whole volume and 
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P (n) = N'-. — (4) 



is the one point distribution function. The function g (ri , r2 ) is the pair correlation function 
(PCF). For homogenous systems (as will be considered here) p (rj ) gives the electrons 



number density p (r j ) = p = N/V {V = system volume) and the PDF and PCF depend 
only on the particles separation: p (ri2) = p^g (ri2) . 

For the homogeneous electron gas in three and lower dimensions, the PDF as well its 
Fourier transform, the static structure factor S (k), have been studied by several authors 
using diverse analytical techniques[6] and also simulations methods. In 3D, we mention ran- 
dom phase approximation (RPA) calculations |10j. diagrammatic ladder approximations [TT] . 
the local-field based method of Singwi, Tosi, Land and Sjolander (STLS)[T2], calculations 
with quantum hypernetted chain equations (QHNC)[13] and also techniques that use pa- 
rameterized PDF, the parameters being determined from known independent theoretical or 
simulation results [Tl]-[l7j. From the side of Monte Carlo quantum simulations in both -the 
variational and diffusion- versions, the work by Ortiz and Ballone[T8] extends in several ways 
previous results of Ceperley and Alder [2], [3]. 

Most of these approaches lean on the variational principle according to which the energy 
of the ground state Eq is a lower bound for the Hamiltonian mean value as calculated using 
any trial wave function ipx {ri ,r2 , ■ ■ ■ , r^r ): 

As trial function a factorized form 

?/^T (ri , r2 , ■ ■ ■ , r^v ) = ^ (ri , r2 , ■ ■ ■ , r^v ) ^0 (ri , r2 , ■ ■ ■ , Fat ) (6) 

is frequently used. Here ipo {ri ,r2 , ■ ■ ■ , rjv ) denotes the system wave function when the 
interactions are turned off (f (rj_,) =0). It is an antisymmetric function under particles 
permutations. We can write 



^o(ri ,r2 ,r^) = ^ (-1)^ P {0i (ri ) , 02 (r2 ) , ■ " " An {rN )} = det[(f)aA^j )] (7) 

p 

where P is the permutation operator that interchanges the particle positions, 0q,- (r^ ) is the 
wave function of an isolated particle and det [0a. {rj )] means the Slater determinant. The 
symmetric factor F (ri , r2 , ■ ■ ■ , r^r ) accounts for the correlations among the particles when 
the interactions are turned on. The A^-body correlation factor can be factorized according 
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to Jastrow[T9]: 

In this work we consider the evaluation of the PDF for the homogeneous 3D electron gas 
starting from a parameterized trial wave function of the form given by Eqs. ([6||8| with the 
parameters obtained by optimizing the system energy by means of a genetic algorithm. 

Genetic algorithms |20j- [2T] form part, together with evolutionary programming [22] . |23] . 
game-playing strategies [21], genetic programming [25] and other related techniques, of a rel- 
atively new class of optimization algorithms which are based on the Darwinian evolution 
principle [26]. In particular, genetic algorithms tackle even complex problems with surprising 
efficiency and robustness. In Physics they have been used in calculations that involve from 
simple quantum systems [27] to astrophysical systems [2S], running through lattice models 
for spin glasses[29], molecules [30] and clusters[31j. More recently[32] we have developed a 
genetic algorithm for the PDF of the one-dimensional electron gas in what, at our knowl- 
edge, is the first application of this kind of algorithm to describe many-body systems in the 
thermodynamic limit.. 

In general, a genetic algorithm is based on three main statements: 

a) It is a process that works at the chromosomic level. Each individual is codified as a 
set of chromosomes. 

b) The process follows the Darwinian theory of evolution, say, the survival and reproduc- 
tion of those individuals that best adapt in a changing environment. 

c) The evolutionary process takes place at the reproduction stage. It is in this stage when 
mutation and crossover occurs. As a result, the progeny chromosomes can differ from their 
parents ones. 

Starting from a guess initial population, a genetic algorithm basically generates consecu- 
tive generations (offprints). These are formed by a set of chromosomes, or character (genes) 
chains, which represent possible solutions to the problem under consideration. At each algo- 
rithm step, a fitness function is applied to the whole set of chromosomes of the corresponding 
generation in order to check the goodness of the codified solution. Then, according to their 
fitting capacity, couples of chromosomes, to which the crossover operator will be applied, are 
chosen. Also, at each step, a mutation operator is applied to a number of randomly chosen 
chromosomes. 
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The two most commonly used methods to randomly select the chromosomes are: 

i) The roulette wheel algorithm. It consists in building a roulette, so that to each chro- 
mosome corresponds a circular sector proportional to its fitness. 

ii) The tournament method. After shuffling the population, their chromosomes are made 
to compete among them in groups of a given size (generally in pairs). The winners will be 
those chromosomes with highest fitness. If we consider a binary tournament, say the com- 
petition is between pairs, the population must be shuffled twice. This technique guarantees 
copies of the best individual among the parents of the next generation. 

After this selection, we proceed with the sexual reproduction or crossing of the chosen 
individuals. In this stage, the survivors exchange chromosomic material and the result- 
ing chromosomes will codify the individuals of the next generation. The forms of sexual 
reproduction most commonly used are: 

i) With one crossing point. This point is randomly chosen on the chain length, and all 
the chain portion between the crossing point and the chain end is exchanged. 

ii) With two crossing points. The portion to be exchanged is in between two randomly 
chosen points. 

For the algorithm implementation, the crossover normally has an assigned percentage that 
determines the frequency of its occurrence. This means that not all of the chromosomes will 
exchange material but some of them will pass intact to the next generation. As a matter 
of fact, there is a technique, named elitism, in which the fittest individual along several 
generations does not cross with any of the other ones and keeps intact until an individual 
fitter than itself appears. 

Besides the selection and crossover, there is another operation, mutation, that produces a 
change in one of the characters or genes of a randomly chosen chromosome. This operation 
allows to introduce new chromosomic material into the population. As for the crossover, 
the mutation is handled as a percentage that determines its occurrence frequency. This 
percentage is, generally, not greater than 5%, quite below the crossover percentage. 

Once the selected chromosomes have been crossed and muted, we need some substitution 
method. Namely, we must choose, among those individuals, which ones will be substituted 
for the new progeny. Two main substitution ways are usually considered. In one of them, all 
modified parents are substituted for the generated new individuals. In this way an individual 
does never coexist with its parents. In the other one, only the worse fitted individuals of the 
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whole population are substituted, thus allowing the coexistence among parents and progeny. 

Since the answer to our problem is almost always unknown, we must establish some 
criterion to stop the algorithm. We can mention two such criteria: i) the algorithm is run 
along a maximum number of generations; ii) the algorithm is ended when the population 
stabilization has been reached, i.e. when all, or most of, the individuals have the same 
fitness. 

In Section III we will apply these ideas to determine the parameters appearing in the 
expression for the 3D electron gas PDF that we propose in Section II. 



II. PARAMETERIZED PDF 



We assume that the trial wave function for the system of A^-spinless electrons with 
Hamiltonian given by Eqs.(T][2) has the form of Eq.(|6]) where we use for the ideal part 
ipo i^i ,r2 , ■ ■ ■ ,rN ) a parameterized generalization of an expression given by Lado[33] and 
for the Jastrow correlation factor also a parameterized expression containing a random phase 
approximation (RPA) pseudopotential[2] . [31] . Specifically we propose 



?/'T (ri , r2 , ■ ■ ■ , r^v ) = exp I [awo (nj) - ^Urpa (nj)] | 



where a and (3 are the parameters to adjust. 

In Eq.(|9]) the ideal gas effective potential Wq (r) is defined 



(9) 



Wo (r) = In go (r) - 



dke~ 



(10) 



Here go (r) and 5*0 (k) are the ideal PCF and structure factor, respectively, whose expressions 
are ^35]. [36]: 



and 



9o [r) 



sin (kpr) — kpr cos (kpr) 
{kprf 



(11) 
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(12) 



with kp the Fermi momentum kp = {3tt p*)^ where p* = pao (ag is the Bohr radius). 
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The RPA pseudopotential, on the other hand, reads 



2p*URPA (k) 



Soik) 



+ 



[So {ky 



Amp*v {k) 

n?k^ 



1/2 



(13) 



where v [k) is the Fourier transform of the interparticle potential v (r). 

It is worth mentioning that in the genetic algorithm me have develop for the ID electron 
gas in Ref.[32], instead of using the RPA pseudopotential in Eq.(|9| we assume that u{r) 
is an unknown function and the algorithm is designed to completely obtain it. Here, a 
parameterized form for the pseudopotential is proposed a priori and the algorithm looks for 
the optimal parameters. 

In principle, to calculate PDF from Eq.([3]) we have to integrate the square of the wave 
function over 3A^ coordinates. To avoid this formidable task use is done of a modified form 
of the hypernetted chain approximation (HNC) for which the PDF is written as a single 
integral equation involving just a pair of particles. We write, ignoring all the elemental 
diagrams[37j. 



g (r) = exp [awo (r) - I3urpa (r) + (r)] 



(14) 



N{r) = p [g (r') -l-N (r')] [g {\v - v' \) - I] dv 



(15) 



where N {r) denotes the sum of nodal diagrams. Eq. (14) without the term with a in the 



exponential has the form of the PCF for a system of bosons (see v.g. Eq.(40) in Ref. 
where the nodal diagrams are denoted D and the elementals diagrams E must be taken 
zero). The term with a adds the ideal part that contains the proper symmetry for fermions. 

Finally, in the variational approach which is implicit in Eq.([5]), we need an expression for 
the Hamiltonian mean value. Making use of Jackson-Feenberg identitv[38] we obtain 



E 

N 



P 



drg (r) 



-— VMn/2 (r)+v(r) 
Am ^ ^ 



(16) 



III. PARAMETERS OPTIMIZATION 



Our problem is to solve for g (r) the integral equation given by Eqs.(14 and 15) with the 
parameters a and /3 determined by demanding that the energy functional E = E {g (r)} 
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given by Eq.(16) be minimum. To this end we use a genetic algorithm. 

We proceed by first generating the initial population, which is formed by Np random 
replicas of the two numbers string (that represents one population individual) ^("),^(^) 
where 7*^"-' G [0, 1] is a random real number (rounded to an established number n of decimals) 
which is assigned to the parameter a. Given a string the encoding consists in 



replacing the sequence of real numbers by a single natural number obtained by putting their 
lecimals ] 

'y2^\..'yi^^ then we have the chain : 



decimals parts one next to the other. Thus if 'y^'^'^ = 0.7!"^ 72"^..7i°'' and 7*^^^^ = 0.7p^ 



(a) (a) (a) (/?) (/?) (B) 

7i 72 •••7; ^7i 72 •••7;'^^ 



and the population is the set 



{( 



In genetic terms, the encoding produces the chromosomic structure of the individuals 
(replica string). The inverse process is called decoding and returns the parameters a and 
P corresponding to each individual. In the decoding we allow the returned parameters be 
multiplied by a constant factor rj > 1 in order that the parameters can take values in an 
interval wider than [0, 1]. We define the fitness of the rth individual as fr = e"^'' where Er 



is the energy calculated using Eq.(16) when g (r) is calculated from Eqs.(14 and 15) for the 
parameters a and (3 obtained by decoding the chromosome structure of the rth individual 
of the population. A solution is reached when f^^l for some individual r in some of the 
successive populations obtained in the evolution process. 

The calculation proceeds by dividing the population of Np replicas into Np/2 couples. 
The couples are randomly chosen by using the roulette wheel algorithm [2Uj. This is done 
by defining the sums F = J2^Zifr and Ss = J2t=i fr = 1, 2, A^p). Then, a random 
number k G [0,F] is generated and the unique index 6 such that Sg-i < k < Ss is picked 
up. 

Once the first generation of replicas (parents) has been generated and divided into cou- 
ples, the second generation (offspring) can be generated by applying the crossover operator 
between the members of each one. At times , some of the members of the new replicas 
generation can be changed by applying the mutation operator. 
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TABLE I: Parameters a and /3 





a 


/3 


1 


2.21TJ 


1.6456 


2 


7.8707 


14.865 


3 


9.6789 


30.526 


6 


14.710 


43.030 


10 


44.974 


18.464 


50 


-6834.17 


128.71 



Given a couple of replicas, the crossover operator is defined by generating a new random 
number c G [0, 1] which is compared with a pre-established crossover probability "p G [0, 1]. 
If c < p, the crossover operator acts by interchanging all the digits from the sth position to 
the end of the replica between the members of the couple. Here s is a random integer such 
that 1 < s < 2n. for example if the couple is 

153280472.. .337 

768325399.. .069 

and s = 4, the new offspring couple will be 
153225399. ..069 
768380472.. .337. 

To apply the mutation operator we first randomly select those offsprings that will mu- 
tate. Then, for each of these offsprings, a gene (a digit) randomly chosen is changed by a 
random integer number I G [0, 9] .The algorithm is stopped when a solution is reached for 
the parameters a and /3. 



IV. RESULTS 



As it is easily seen, the electron gas is completely determined by giving just its density p 
or, as it is custom in many-body theory, the Wigner-Seitz radius defined = [3/47rp]^''^. 
Here we have applied the procedure described above, to an electron gas at metallic densities: 
1 < '"s < 10 and also at = 50. We use in our calculations = 150 and n = 5. The factor 
f], in turn, is moved in each case to give reasonable values for the parameters a and /3. In 
Table [T] we show the parameters a and (3 obtained for the diverse values of considered. 

9 



U,U-| 1 1 1 1 1 1 1 1 

12 3 4 

r 

FIG. 1: The functions g (r) for the homogenous electron gas at rs=l,2,3,6,10 and 50 obtained in 
this work. 

Note the change in the parameters tendency at = 50. By the way, we were unable to 
reach convergence for values of greater than this one. 

Figs. 1 to 3 show the corresponding g (r)'s. In Fig. 1 we put all together the curves 
we have obtained . We observe the characteristic features of the electron gas: when the 
density increases (r^ decreases) the behavior tends to that of an ideal paramagnetic gas of 
fermions with contact value 1/2 and rapidly going to the asymptotic value 1. When the 
density decreases the correlation functions become more structured showing, in particular, 
a more pronounced Coulomb hole near contact. 

In Figs. 2 and 3 comparison is done of our results with those obtained from variational 
as well as diffusion Monte Carlo simulations performed by other authors p!H] .113]. . A first 
remark is the existing differences between variational and diffusion results particularly for 
low values of r^. Also must be noticed the good agreement of the results of this work with 
those obtained from diffusion Monte Carlo calculations. 

It is worth mentioning that the time to obtain one of our curves by running the complete 
algorithm with a Pentium IV is tipically of the order of 50 hours for the metallic densities. 
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FIG. 2: The functions g{r) for the homogenous electron obtained in this work compared with 
variational and diffusion Monte Carlo results, (a): r5=l; (b): rs=3. 
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FIG. 3: The functions g{r) for the homogenous electron obtained in this work compared with 
variational and diffusion Monte Carlo results, (a): rs=10; (b): rs=50. 
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